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We recently derived a very accurate and fast new algorithm for numerically inverting the Laplace 
transforms needed to obtain gluon distributions from the proton structure function _F,7 p (:r, Q 2 ). We 
numerically inverted the function g(s), s being the variable in Laplace space, to G(v), where v 
is the variable in ordinary space. We have since discovered that the algorithm does not work if 
g(s) — ¥ less rapidly than 1/s as s — > oo, e.g., as 1/s 13 for < /3 < 1. In this note, we derive 
a new numerical algorithm for such cases, which holds for all positive and non-integer negative 
values of /3. The new algorithm is exact if the original function G(v) is given by the product of a 
power w' 3-1 and a polynomial in v. We test the algorithm numerically for very small positive /3, 
/3 — 1CP 6 obtaining numerical results that imitate the Dirac delta function 5(v). We also devolve the 
published MSTW2008LO gluon distribution at virtuality Q 2 — 5 GeV 2 down to the lower virtuality 
Q 2 — 1.69 GeV 2 . For devolution, (3 is negative, giving rise to inverse Laplace transforms that are 
distributions and not proper functions. This requires us to introduce the concept of Hadamard 
Finite Part integrals, which we discuss in detail. 



I. INTRODUCTION 



In an earlier note [T] we developed an algorithm to numerically invert Laplace transforms in order to find an analytic 
solution for gluon distributions, using a global parameterization of the proton structure function, F2 P (x,Q 2 ) and a 
LO (leading-order) evolution equation for F^" 9 . However, when we went to NLO (next-to-leading order) in the strong 
coupling constant a s (Q 2 ), we have discovered the algorithm failed badly. Detailed investigation showed that the 
cause of the problem was that this <7(s)-the Laplace transform of our desired NLO gluon distribution G(v), where 
v = ln(l/a;) — went to less rapidly than 1/s as s — > oo, where s is the Laplace space variable. The purpose of this 
note is to derive a new and exact algorithm for such cases, which can be modeled by Laplace transforms of the type 

1 M ~ x h 

fe=0 

for all values of positive (3 and all values of non-integer negative /?. We note that for f3 — 0, —1, —2, . . ., the inverse 
Laplace transforms are the distributions S(v),6'(v),6"(v), . . ., the Dirac delta function and its derivatives, which are 
not true functions but rather are distributions. 

In Sections [TTJ |III| and |IV| we derive exact numerical Laplace inversions for originals of the form 

A/-1 , 

but now generalized for all positive values of f3 together with all negative — but non- integral — values of f3. In this 
context, exact means calculation to arbitrary numerical precision, given a symbolic program such as Mathematica [2 
which also calculates numerically to arbitrary accuracy. If the function G(v) is well-approximated by G(v), one can 
evaluate G(v) to arbitrary accuracy. 

In Section [Vj we show that we can successfully reproduce the equivalent of a numerical delta function, using a very 
tiny positive ft in Eq. M. To illustrate the new method, we will numerically invert g(s) — i/s 1 / 1000000 ^ the Laplace 
transform of v~ 1+1 / 100 ™ 00 /F(l/1000000) (a numerical surrogate for a Dirac 5 function) and test its accuracy. 
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In Section VI we solve a real physical problem, the devolution of the published LO MSTW 2008 [3 gluon distribution 
from the virtuality Q 2 = 5 GeV 2 to Q 2 = 1.69 GeV 2 (the squared mass of the c quark). This calculation involves 
a rather large negative value of /?, and consequently, calculation of a distribution, rather than a function. We must 



numerically compute a convolution integral in Eq. (62), i.e., 



K GG (w)G (v - w) dw, 



(3) 



where the kernel K GG (w) is given by Eq. ^ with negative ~ —0.5, i.e., K GG (w) is a distribution (not a function 
in the usual sense) in w about half way between 5(w) and S'(w)), and thus the above integral is divergent. To obtain 
a convergent result, we must replace the Riemann integral sign J Q U in Eq. (|3j) by the Hadamard Finite Part integral 
sign , which introduces the regularization obtained by using the Hademard "parte finie" (Finite Part) integral [4], 
discussed in depth in Appendix [C] 



II. NUMERICAL INVERSION OF LAPLACE TRANSFORMS 



Let g(s) be the Laplace transform of G{v). The Bromwich inversion formula for G(v), which we call the original 
function, is given by 

G(v) = C- l W)'A = ^ / dsg(s)e vs , (4) 

where c is a real constant such that the integration contour lies to the right of all singularities of g(s). We will assume 
that we have made an appropriate coordinate translation in s space so that those singularities all lie in the left-half 
complex plane, and take c = 0. Our goal is to numerically solve Eq. The inverse Laplace transform is essentially 
determined by the behavior of g(s) near its singularities, and thus is an ill-conditioned or ill-posed numerical problem. 

In this note we present a new algorithm that takes advantage of very fast, arbitrarily high precision complex 
number arithmetic that is possible today in programs like Mathematica [5] , making the inversion problem numerically 
tractable. 

First, we introduce a new complex variable z = vs and rewrite Eq. Q as 

We assume that the form of the Laplace transform g[s) can be approximated as 

1 M_1 b 

g(s) « g(s) = J2 ( 6 ) 

fc=0 

corresponding to an original function G(v) which can be approximated as in Eq. ^ as a factor v^~ x times a polynomial 
of order M — 1 in v, i.e., 

M-i , 



The sum contains M coefficients bj. 



In our earlier paper jT], we proceeded to make a rational approximation for the exponential e z , under the tacit 
assumption that g(s) went to sufficiently rapidly as s went to oo, and used this to evaluate the integral in Eq. jjjj. 
This replacement is not useful numerically if g(s) goes to too slowly for s — ¥ ±ioo, e.g., as with /? < 1 [5]. To 
generalize for all possible /?, we now rewrite Eq. ^ as 



+i oo 

(8) 



G(v) = J- / dzg (-) z?- 1 



e 



Our new algorithm uses a different treatment of the factor in square brackets. 
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A. The case of P = 1 

In order to understand the more complicated case of f3 ^= 1, we briefly review here the procedure followed for the 
solution (3 — 1, detailed in Ref.p]. When j3 = 1, Eq. ([6| is given by 

A/-1 



*(•) - E I w 

and Eq. Q is given by 

G(v) = ^r.>^Q[e1. (10) 



s f — ' s" 

k=0 



+i oo 



In the earlier paper, we approximated the square bracket [e z ] by a rational function defined as a ratio of polynomials 
with the numerator a polynomial of order 2N — 1 and the denominator a polynomial of order 2N . This gives 



2N 



E 



z - a,- 



(11) 



Inserting this into Eq. using /? = f , we obtained 



+ 2 OO 



2 AT 



The key observation for numerical purposes was that, by closing the contour of integration in the right half of the 



complex plane, possible with the expression in Eq. \ll\ but not possible for e z itself, the integral in Eq. ( 12 ) could be 
evaluated simply as a sum of the residues of the integrand at the poles otj. No numerical integration was necessary, 
and there were no contributions from singularities of g{z/v) which lie entirely in the left half plane. 

As we showed in Ref. [I], the requirement that the result be exact for all inverse powers l/s" with 1 < n < AN, i.e., 
that M = AN in Eq. Q, was sufficient to determine the AN complex parameters on, Wj uniquely, and furthermore, to 



show that the expression in Eq. (II ) was equal to the Pade approximant P(e z ; 2N — 1, 2N) of e z . This is defined as 
the ratio of polynomials in z of orders 2N — 1 and 2N which, when expanded, exactly reproduces the first AN terms 
in the Maclaurin expansion of e z . Conversely, the use of the (2N — 1, 2N) Pade approximation to e z automatically 
gave exact results for powers of l/s in the stated range, or powers of v in the original function up to v ^. 

The Pade approximant converges to e z to arbitrary accuracy for N sufficiently large. Its use instead of e z is justified 
for a given N if the error in the approximation to e z is sufficiently small, and g(z/v) vanishes sufficiently rapidly for 
s — > ±ioo, that the approximation is valid over the region in s that contributes significantly to the integral. 



It can be shown from the properties of the Pade approximant that the poles ctj and weights u>j in Eq. ( 11 ) have the 
following properties: 

1. The poles a,j (the zeros of the Pade denominator) are all distinct and appear in complex conjugate pairs which 
we label as (aj,aj + i), j odd, with <Xj+i = ctj the complex conjugate of aj. The poles have Re a } > for all j, 
so are all in the right-hand half of the complex plane. 

2. The weights u)j also appear in complex conjugate pairs so there are N distinct pairs of the complex numbers 
(bjj,a,j), such that the sum of the fc th pair is real, 



*3+l 



S Real, j = odd. (13) 



3. The integrand vanishes faster than l/R as R — > oo on the semi-circle of radius R that encloses the right portion 
of the complex plane, since the approximation vanishes as l/R and g(z/v) also vanishes for R — > oo. 



To evaluate the integral in Eq. ( 12 ) using these properties, we formed a closed contour C by completing the inte- 
gration path with an infinite half circle in the right portion of the complex plane, where g(z/v) has no singularities. 
It is important to note that this contour is a clockwise path around the poles aj of Eq. (12), which arise because 
we replaced the square bracket of Eq. (jllj) by Y^j=i u j/{ z — otj). What we need is the negative of this path, i.e., the 
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contour — C which is counterclockwise, so that the poles are to our left as we traverse the contour — C . Accordingly, 



Eq. (12 1 was rewritten as 



G(v) 



2mv j r^iv) E 

J ~ C 3=1 

2N 
3=1 

^E Re Kf) Wj 



(14) 
(15) 
(16) 



3 = 1 



To obtain Eq. ( 15 ), we used Cauchy's theorem to equate the closed contour integral around the path — C to 2ni times 
the sum of the (complex) residues of the poles. Since the contour —C restricts us to be on the right of any singularities 
of g(s), no poles of g(z/v) are enclosed; only the 2N poles on are inside the contour. To obtain our final result for 
G(v) in Eq. ( 16 ), we used the properties cited above of the complex conjugate pairs in Eq. ( 13 ). Taking only their real 



part and multiplying by 2, we have simultaneously insured that G(v) is real, yet only have had to sum over half of 
the residues. 



B. Generalization to /3 7^ 1 



The situation is more complicated for < /3 < 1. The Laplace transform g(s) is then of the form in Eq. 



M-l 



5( s ) = i E 



b k 



k=0 



(17) 



with /3 non-integer, and the original algorithm in [T] fails numerically. 
To handle this case, we first rewrite the inversion formula in Eq. Q as 



G(v) 



2iriv 



+i 00 



dzg ( — I 



P-i 



rfi-l 



(18) 

and causes 
there is no 



The function ■/P~ 1 g(z) is of the form in Eq. ffify, with the original non-integer value of /3 replaced by 1, 
no difficulty. The problem arises from the term in square brackets, [•] = e z /z@~ 1 . For (3 non-integer, 
Maclaurin series about z = 0, hence no Pade approximant or rational approximation for this factor. 

It is still the case that the previous method works for any positive integer value of j3, j3 = n > 1, for 2N sufficiently 
large, a result that follows from the replacement of the first n — 1 coefficients bk in Eq. ([9| by zero. This suggests that 
the replacement of [•] by a rational function may still be useful for (3 non-integer. We therefore replace the bracket on 
the right-hand side of Eq. ( 18 1 by 



yj-l 



2N 

E 



P~l( V 



(19) 



a rational function of z which retains the form which works for integer values of j3. We stress that we do not regard 
this replacement as giving an adequate approximation to the term in brackets; the behavior of the two functions is 
quite different for z — > for non-integer /?. 

Instead, we change our emphasis from approximation to exactness, and require that the coefficients aj and ujj be 
chosen such that the new expression gives exact results when integrated with all inverse powers 1/ s n with 1 < n < AN . 

In the case of /3 = 1, the condition of exactness required that the rational function be the Pade approximant of 
e z ; conversely, the use of the Pade approximation for e z led automatically to exactness. We will show here that an 
appropriate choice of the coefficients a,- and ujj is possible for non-integer /3, and that the function on the right-hand 



side of Eq. (|19| is the (2iV-l,2iV) Pade approximant for the function p(z) = iFi(l,f3,z)/T(f3) where i-Fi(l, f3, z)/T(f3) 
is the Kummcr confluent hypergeometric function. 

This result allows us to readily obtain exact inverse Laplace transforms for functions with g(z) of the form in Eq. ([T]) 
or Eq. (17), corresponding to original functions that can be approximated by finite series of the form in Eq. ([2]), thus 
justifying the seemingly arbitrary replacement a posteriori. 
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With this replacement, Eq. ( 18 1 can be evaluated as was done before, by closing the integration contour in the right 



half complex plane and using the Cauchy residue theorem, giving 

r+iao 2JV 

2iriv 

2N 



G(v) 



, 2N 



(20) 
(21) 



III. CANONICAL EQUATIONS 

The task now is to find the appropriate AN poles otj and weights uij, j — 1,2,..., 2N, such that the expression 
in Eq. (21 1 gives the exact inverse in Eq. (21 1 for all functions g(z) of the form in Eqs. |l]) or (IT). This requires that 



4N-1 



G(v) 



^ E 

k=0 
4N-1 2N 



T(p + k) 



-, 2N 

v z — ' V V J 



3 = 1 



= - J2 vP+k ~ lb -~ : 



k=0 3 = 1 



(22) 
(23) 



Equating the coefficients of the arbitrary parameters bk, we find that the ay and uij must satisfy the set of AN 
simultaneous equations 



2N 



r(/? + fc)]T 



3 = 1 a 3 



-1, fc = Q,l,. 



,47V- 1. 



(24) 



As is necessary, these canonical equations are independent of v. The solutions of the AN equations determine the 
2N weights and 2N poles. However, the canonical equations are ill-posed for numerical purposes, so that direct 



solution of Eqs. (24|) for weights and poles is basically impractical for even modest 2N. Fortunately, the poles and 
weights are readi 



y found by other means, as will be shown in the following Section [TV 



IV. DETERMINATION OF THE POLES ay AND THE WEIGHTS ojj FOR ARBITRARY /3 
A. Connection with Pade approximates for confluent hypergeometric functions 



We will now rederive the canonical equations of Section |III| by a completely different technique which makes it a 
simple numerical task to calculate accurately and quickly the 2N poles aj and 2N weights uij that satisfy the AN 



canonical equations of Eq. ( 24 ) . 



Let us now consider a generic term of Eq. (17), 

9k(s) = - 



1 



0,1, 



,M-1. 



The exact inverse Laplace transform of §k( s ) is given by 

1 



c- 1 



cP+k 



:v = l 



r(/3 + k) ■ 



(25) 



(26) 



where, for simplicity, we have evaluated it at v = 1. 

We next revisit the replacement of [e z /z^ 1 ] by a rational function which we made in Eq. (20), defining the rational 
function as V(z), 



2 A 



V{z) 



E 



(27) 
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This can obviously be written as the ratio of polynomials in z of orders 27V — 1 and 27V containing a total of 47V 
coefficients, 27V in the numerator and 27V in the denominator when the first term in the denominator is fixed to equal 
1; this is just the form of a (27V — 1, 2N) Pade approximant. 
We now rewrite Eq. (20) in terms of V{z) as 



1 



2iriv 



:g(l)z^Viz) 



k\ dz k ^ ^ 



(28) 
(29) 



where in Eq. (29) we have substituted the Maclaurin expansion of the rational function V(z). 



We next examine the conditions imposed on Viz ) by the requirement that the result of the int egra tion in Eq. ( 28 ) 
be exact for functions §(s) of the form in Eq. (25). When we replace g(z/v) by l/z l3+k in Eq. (29)_for v = 1 and 



close the integration contour to the left rather than the right, the result of the integration of Eq. ( 29 ) is the residue 
of V{z)/z k+1 at z = 0, namely the coefficient of z k in the Maclaurin expansion of Viz) around 2 = 0. The condition 
that the result agree with the exact evaluation of the integral in Eq. ( 26 1 thus requires that 



1 d k V 
k\ dz k 



(0) 



r(/3 + k) 



0,1,2, 



AN -I. 



(30) 



We can now adjust the 47V coefficients in Viz) so that these conditions are satisfied for the first 47V terms in the 
Maclaurin expansion of Viz). 



Let us now extend the sum on the right hand side (r.h.s.) of Eq. (29) to infinite N , and define a new function 



piz) = 



r([3) r(/? + i) r(p + 2) res + 3) 



E 

k=0 
1 



Tif3 + k) 



(31) 
(32) 

(33) 



where /3, z) is the Kummer confluent hypergeometric function. 

It is readily seen that V {z) is the Pade approximant of order (27V — 1,27V) for the function p(z) = ii<i(l,/3, z)/T(0): 
it is the ratio of polynomials of orders 2N — 1 and 27V which, when expanded, reproduces the first 47V terms in the 
Maclaurin series for p(z). This is precisely the definition of the Pade approximant for pjz), so Piz) = P(p(z), 27V — 
1,27V). The poles ctj and weights Wj, j = 1,2, ...,27V, needed for the result in Eq. (21) to be exact for giz) a 
polynomial of degree 47V — 1 in 1/z, or G(v) a polynomial of degree 47V — 1 in v, are given by the 27V zeroes of the 
denominator and the 27V residues of the Pade approximant at the poles. 

Although we are not aware of a formal proof except for the case [3 = 1, the poles ctj are found in practice to appear 
only in complex conjugate pairs in the right half of the complex plane, The residues also appear in conjugate pairs. 
Using this information, we rewrite Eq. (21 ) as 



G(v) 



- 2 x>K^) 

v r — ' l \ v / 



(34) 



where we sum only over the pol es in the upper half plane. 

We emphasize again that Eq. (34) becomes an exact statement if G(w) is the product of v^~ x times a polynomial of 
order 47V — 1. Obviously, if an original function G(i>) can be adequately approximated as the product of v^^ 1 times 
a polynomial of order 47V — 1, we can then approximate G(v) by G(v) and write 



GM 



(35) 



i=l 



We have found the algorithm to work very well in practice, even for fairly small values of 27V. 
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Revisiting Eq. ( 19 1, we see that we have completely justified the replacement of the bracket by the sum for arbitrary 
(3 > 0, i.e., 

27V 

■+E rz? r » (36) 



0-1 



(z — Oj) 



not as an approximation for [■], but as a method to obtain exact or essentially exact numerical results for a large class 
of functions. 

We can readily calculate the 27V pole positions aj (the zeroes of the denominator of P(p(z),2N — 1,27V)), and 
the 27V weights uij (the residues of P(p(z), 27V — 1, 27V) at the poles) to arbitrary accuracy using a program such as 
Mathematical [2], which can use arbitrary-accuracy complex arithmetic. 

We will show in Appendix |A| that the Pade approximant P(p(z),2N — 1,27V) can actually be obtained in closed 
form and calculated rapidly, without the necessity of calculating 27V — 1-fold derivatives of the function p(z). As we 
already pointed out, since both the cu and the ujj occur in complex conjugate pairs, we only have to calculate half 'of 
them to use the relationship in Eq. ( 34 ) . 

A concise inversion algorithm in Mathematica that rapidly and accurately implements Eq. ( 16 ) is given in Appendix 

m 

Finally, we note as an aside that the function p(z) = i-Fi(l, /?, z)/T(/3) is particularly simple for some special cases. 
For /3 = 1, it becomes 

(37) 

z 1 / 2 e z , the 



p(z) = e\ [3= 1, 



II A 



and the construction above reproduces the results found in Sec 
example used for illustration of a different algorithm in [5], p(z) becomes 

1 

^pK Res— >-oo 

-* 2 alt. 



p{z) 



z 1 l 2 e z eii{^) 



and [J. For (3 = 1/2 and e z / z^ 1 
■mes 

zVV, = 1/2, 



(38) 



where erf(z) is the error function, given by erf (z) = J Q e 

More generally, the leading term in the asymptotic expansion of p(z) for Rez — > oo is just e^/z^ -1 , the factor in 
square brackets in Eq. (18 1. However, there are additional asymptotic terms which increase less rapidly for Rez — > oo. 
As already remarked, the behaviors of p(z) and e z /z^~ x are quite different for z — > 0; they also differ for z — > ±oo. 



B. Analogy to Gaussian integration routines 



We now slightly rewrite Eq. ( 34 ) as 



1 r +io ° { e z \ J-L 
VG{v)= 2^ij_- F{Z) \~z^) ^«- 2 E Rc [^)^]- ( 39 ) 



where 



and the weights Wi are 



F(z) = z^ 1 ~g[-) (40) 



W i = uj i /a?- 1 . (41) 

The result is exact if F(z) is a polynomial in 1/z of degree up to 47V — 1 because of our adjustment of the 27V poles 
cti (the location of the zeros of the Pade denominator) and the 27V residues atj. 

The factor e z /z^ _1 in the integral of Eq. (39) evidently plays the same role in inverse Laplace transformations as 



the weight function u(z) plays in numerical integration done by Gaussian quadrature, where the definite integral I in 
the interval {— 1,+1} is approximated by 

N 

(42) 

3= ' 

The Zi in Eq. ( [42] ) are the TV zeroes of the appropriate orthogonal polynomial P/v(z) in the interval {—1, +1} for the 
weight function u(z), and the weights wj are the Christoffel numbers of the P^{zj). This result is exact if f{z) is a 
polynomial of order less than or equal to 27V — 1. since there are 27V coefficients, the TV zeroes Zi and the TV weights 
Wi to adjust. 



/ u{z)f{z)dzK,^f{ Zj )w 3 . 

J -l 1=1 



With the interchange of poles and zeros, analogy between Eqs. (35) and (42 ) is clear 



V. A NUMERICAL APPROXIMATION OF A DIRAC DELTA FUNCTION 

The Dirac delta function 5(v) is a distribution defined by its property that for any smooth test function f(w) 

rb 



f(v)S(v)dv = /(0), (43) 

^0 

when b > 0. 

When /3 = in Eq. (JlJ, the result is G(v) — S(v), which of course is impossible to write as a true function. We now 
charge our numerical Laplace inversion algorithm with the daunting task of finding an accurate numerical inversion 
of a Laplace transform of a close approximation of the Dirac delta function S(v). We approximate it by 

6{v) a G(v) = w- 1+1 / 100000 7r(l/1000000), (44) 

which uses the very tiny positive value of (3 = 1/1000000 = 10~ 6 for our delta function approximation. For f(v) = 1, 
J h G(v) dv is given by 0.99999 for the upper limit b = 0.0001 and by 1.00001 for b = 10000, compared to the expected 
value of 1, while for the test function f(v) — sin(i>), the integral sm(v) G(w) dv is given by 1 x 10 -10 for b = 0.0001 
and 1.6 x 10 -6 for b = 10000, compared to 0, thus showing that G(v) is a good numerical approximation to the Dirac 
delta function 5(v) over an enormous range of v. The Laplace transform of G(v) is given by 

g(s) = L[G{v,s] = ^4^. (45) 

To illustrate the accuracy of the numerical inversion routine we give in Appendix[XJ we plot in Fig. [T]thc numerical 
inversion of g(s) = i/g 1 / 1000000 } the Laplace transform of our numerical approximation of a S function. The solid 



curve is the exact function G(v) = i;~ 1+1 / 1000000 /r(l/1000000) from Eq. (44 1, and the (red) dots are a result of using 
the algorithm in Appendix [Aj using Mathematica [2]. As expected, the agreement is exact: in this case, exact means 
a fractional numerical precision of O(10~ 30 ) in the entire v range. 

In practice, in order to use our algorithm one needs to know that the original function, the inverse Laplace transform, 
is well-approximated by a power law times a polynomial, as well as the numerical value of (3. Normally we only 
have a numerical g(s), so even knowing that the original function is well-approximated with a power law times a 
polynomial, the actual numerical value of /3 is not known. However, if we can evaluate g(s) numerically to arbitrary 
precision, we can determine P by taking consecutive closely spaced numerical values of g(s) at exceedingly large values 
of s (thus simulating an asymptotic expansion of g(s)) and then fitting these numerical results to a power law in s. 
Of course, this degrades the numerical accuracy of /3 considerabl y, b ut is quite viable for practical situations. For 



example, by fitting the function a/s^ to numerical output of Eq. (45) using 30 digit accuracy, we obtain a value for 
j3 that is sufficiently accurate to give a relative accuracy of 0(3 x 10~ 17 ) in inversion, more than necessary for most 
numerical applications. What is far more critical is knowing that the inverse transform (the original) is adequately 
approximated by a power law in v multiplied by a finite polynomial in v. 



VI. GLUON DEVOLUTION USING LEADING-ORDER SINGLET DGLAP EQUATIONS 

In this Section we will use the work of Block, Durand, Ha and McKay 7 (BDHM), who derived analytic decoupled 
singlet solutions to the leading-order (LO) Dokshitzer-Gribov-Lipatov-Altarclli-Parisi (DGLAP) evolution equations 
[8TH0j. We will use their gluon solution to solve the physical problem of gluon devolution — a much more difficult 
problem than gluon evolution. For details of these coupled integro-differential equations and their solutions, see Ref. 
[7j. Here we will devolve the published LO MSTW2008LO :3! 2008 gluon distribution at the virtuality Q 2 =5 GeV 2 
down to Q 2 =1.69 GeV 2 , the square of the c quark mass used by those authors. 



A. Decoupled gluon solution to the LO DGLAP equation 

BDHM first rewrite the standard LO singlet DGLAP equations, normally written in terms of the virtuality Q 2 
and the variable Bjorken-x, in terms of new variables v = ln(l/x), w — ln(l/z). After introducing the notation 
F s (v, Q 2 ) = F s (e~ v , Q 2 ), G(v, Q 2 ) = G(e~ v , Q 2 ), they find that the singlet DGLAP equations have now been written 
in a form such that all of the integrals in the LO DGLAP equations are manifestly seen to be convolution integrals. 
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FIG. 1: A plot of the numerical Laplace inversion of g(s) = i/gVioooooo and the exact original f unc tion G(v) = 
u- 1+1/1000000 /r(1000000) vs. v. The solid curve is G{v) = u~ 1+1/1000000 /r(l/1000000) and the (red) dots are the numer- 
ical inversion of i/ s 1 / 1000000 ] US mg the algorithm given in Appendix |a| The fractional error of each point is 0(1CP 30 ). 



Introducing Laplace transforms allows BDHM to factor these convolution integrals, since the Laplace transform of a 
convolution is the product of the Laplace transforms of the factors, i.e., 



C 



F[w]H[v — w] dw; s 



= c 


[I 







F[v — w]H[w] dw; s 



C[F[v];s]xC[H[v];s]. 



Defining the Laplace transforms of F s (v,Q 2 ) and G(v,Q 2 ) in Laplace space s as 



/(s,Q 2 ) ee C[F s (v,Q 2 );s 
g(s,Q 2 ) = C[G(v,Q 2 );s] ■. 



F s (v,Q 2 )e- sv dv, 
G(v,Q 2 )e~ sv dv 



(46) 

(47) 
(48) 



they find two coupled first order differential equations in Q 2 in Laplace space s that have Q 2 -dependent coefficients, 
which are 



df 



dlnQ 2 



(s,Q 2 ) 



a s (Q 2 ) 

47T 



$ f (s)f(s,Q 2 



99 (s,Q 2 ) = ^?% 3 (s) 5 (s,Q 2 ) 



<91nQ 



47T 



a s (Q 2 ) 

47T 
47T 



&f(s)g(s 1 Q 2 ) 1 
Q g ( S )f(s,Q 2 ). 



The s-dependent coefficient functions $ and are given by 



$ / (s) = 4- - 
e f (s) = 2n f 



' + -m + 2 + 1) + 7s) 



3 V s + 1 s + 2 
1 2 2 



s + 1 s + 2 s + 3, 

2 1 



33 -2nt /l 
= — 1 + 12 



1 



8 /2 



s s + 1 s + 2 s + 3 
1 



- - ^(S + 1) - JE 



3 Is s + 1 s + 2/ ' 



(49) 
(50) 

(51) 
(52) 
(53) 
(54) 
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where ip(x) is the digamma function and je = 0.5772156 ... is Euler's constant. Here a s (Q 2 ) is the running strong 
coupling constant, and for MSTW2008LO [3] is given by the LO form 



a s (Q 2 ) = 



-in 



(55) 



(11 -§n/) HQ 2 /A 2 lf )' 

with rif the number of quark flavors. The value of the LO version of a s at Q 2 = 1 GeV 2 was used in MSTW2008LO 
as a fitting parameter. The QCD parameters A3 was adjusted to reproduce the fitted value at Q 2 = 1 GeV 2 ; A4 and 
A 5 were fixed so that a s is continuous across the boundaries at Q 2 — M 2 and M 2 where rif changes at the masses of 
the b and c quarks. 

The solution of the coupled equations in Eq. ( 49 ) and Eq. ( 50 ) in terms of initial values of the functions / and 
g, specified as functions of s at virtuality Q 2 ,, is straightforward. The Q 2 dependence of the solutions is expressed 
entirely through the function 



t{q\qD 



1 

4-7T 



Q 1 



a s {Q' 2 )d\uQ 



12 



(56) 



With the initial conditions /o( s ) = f( s , Qo) an d 5o( s ) = g{ s i Qo), BDHM find the decoupled gluon solution in Laplace 
space s is given by 



g{s,r) = k gg (s,T)g {s) + k gf (s,T)fo{s). 
where the coefficient functions in the solution are 

sinh (§.R(s)) 



kgg(s, 7") 



cosh ( - 



(2^) 



R(s) 



($ f (3)-$ g (3)) 



R(s) 



(57) 

(58) 
(59) 



with R(s) = y($/(s) - <f> g (s)y + 40 f (s)O g (s). BDHM now define two kernels Kqf and Kqq, the inverse Laplace 
transforms of the fc's, i.e., 



K G g{v,t) = L 1 [kgg(s,T);v], K gf (v,t) = C 1 [k gf (s, r); v}. 



(60) 



It is evident from Eqs. (56) and (59 1 that Kqf vanishes for Q 2 — Qq where t(Q 2 ,Qq) — 0. It can also be shown 
without difficulty that for r — 0, Kgg( v , 0) = 6(v) and that Kgf(v, 0) = 0. 

The initial boundary conditions at Qq are given by F s q(x) = F s (x 1 Qq) and Gq(x) = G(x,Qq). In u-space, 
F s q(v) = F S Q(e~ v ) and Gq(v) = Go(e~ v ) are the inverse Laplace transforms of /o(s) and go(s), respectively, i.e., 

F sQ (v) = r 1 [/oW;»] and G (v) = C^lgo^v]. (61) 
Finally, BDHM writes the LO gluon distribution solution in v-space in terms of the convolution integrals as 

G(v,Q 2 ) = / K GG (w,T)G (v-w)dw+ / K GF (w,T)F s0 (v~w)dw. (62) 
Jo Jo 

The LO gluon solution in Bjorken-a; space, G(x,Q 2 ), is then easily found by going from w-space to x space via the 
transformation, G(x,Q 2 ) = G(\n(l/x),Q 2 ). 



B. Devolution of the LO MSTW2008LO gluon distribution G(x,Q 2 ) from Q 2 = 5 GeV 2 to 1.69 GeV 2 



In order to use Eq. (J62J) to devolve from the leading-order MSTW2008LO gluon distribution [3] at Q 2 = 5 GeV 2 
to Q 2 — 1.69 GeV 2 — the MSTW value at the square of the c quark mass — we must calculate the singular kernel 
Kgg(w,t) at a negative value of r, i.e., r = —0.0332005. Asymptotically expanding k gg (s,T) in s (see Eq. (58l), we 
can write it as 



kgg( s > T ) 



-12t (33-2n/)/3r 
e 1 



(63) 
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which shows that as s — > oo, k gg (s,t) -> s-P, with P = 12r = -0.398406. Since j3 is negative, the kernel K gg (w,t) 



in the convolution integral L K GG (w, t)G {v — w)dw in Eq. (62) is a distribution, and not a function. Hence, we 

must replace this integral with the Hadamard Finite Part integral, K GG (w, t)Gq(v — w) dw, derived in Appendix 
[C| and discussed fully there. 

Using the Mathematica algorithm developed in Appendix |Aj with g= fc S9 (s,r) from Eq. (58), twoN=20, beta=- 
0.398406, precision— =100, the numerical inverse Laplace transform of the kernel K GG (w, t) that we obtained was 
adequately least squares fitted — using the model of Eq. ^ — by 

8—1 32 
w p \ - 



K gg (w,t) = — ^B z w\ (64) 

thus determining the 33 coefficients Bi. Next, we wrote the devolved gluon distribution as 

K GG (w,T)G (v-w)dw+ / K GF {w,T)F s0 (v-w)dw, (65) 



evaluating the Hadamard Finite Part integral (the first integral in Eq. (65), involving the kernel K GG ) using the 
Mathematica algorithm in Appendix [C] while doing a straightforward evaluation of the second (Reimann) integral, 
using /3 = 1 in our Laplace inversion algorithm. Finally, we returned to Bjorken-x space, computing G(x, Q 2 = 1.69) 
using the substitution x = e~~ v . 

Shown in Fig. [2] as the (red) dots are our numerical inversion devolution results for G(x) at Q 2 — 1.69 GeV 2 , 
devolved from Q 2 = 5 GeV 2 , compared to the published MSTW2008LO gluon distribution at Q 2 =1.69 GeV 2 [3]. The 
numerical agreement is excellent, with the fractional error at the smallest x values being 0(1 x 10 -4 ), the accuracy 
with which the gluon distribution is given on the Durham web site [TT]. When we devolve from Q 2 = 10 GeV 2 , the 
fractional error becomes 0(6 x 10~ 4 ); devolving from Q 2 = 15 GeV 2 , the fractional error is 0(3 x 10 -3 ); evolving 
from Q 2 = 20 GeV 2 , the fractional error degrades to 0(2 x 10 -2 ). The decrease in accuracy with the greater range 
of devolution is related to the uncertainties in the initial distributions, which are given only numerically on a grid in 
x. The resulting uncertainties in fits to those distributions grow essentially exponentially under devolution. 



VII. CONCLUSIONS 



By numerically evaluating the original G(v) from the Laplace transform g(s) as 

N 



G(«) 



i=l 



(66) 



using the short Mathematica function given in Appendix |Aj we have achieved an exact numerical solution for inverse 
Laplace transforms whose originals are of the form 



G(v) 



„B-1 A£-l 



B lV l , M < AN. 



(67) 



When G(v) is adequately approximated by the r.h.s. of Eq. (67), we obtain an excellent numerical approximation. 

We note in passing that when (8 = 1, this algorithm also replaces our earlier numerical Algorithm I pQ with a 
slightly more efficient one. 

As an example of a very difficult numerical problem, we used the algorithm of Appendix [X] to invert numerically 
a Laplace transform of an original function which is an excellent approximation to a Dirac 5 function. We show in 
Fig. |lj the numerical inverse Laplace transform of g- 1 / 1000000 together with the exact answer, to demonstrate the 
algorithm's inherent accuracy. A fractional accuracy of O(10 -30 ) was achieved. 

For our second example — a real physical problem — we accurately devolved the published MSTW2008LO gluon 
distribution [3] from a virtuality Q 2 = 5 GeV 2 down to Q 2 — 1.69 GeV 2 , achieving a fractional accuracy of O(10~ 4 ). 

Although our inversion routine was originally developed to calculate the inverse Laplace transforms needed in work 
on the evolution of gluon distributions, it quite general and has a wide variety of potential applications , e.g., in the 
solution of both integral and differential equations. 
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Dots: Devolution of MSTW LO G from 5 to 1 .69 GeV^ 



25 - 



20 - 



x 
O 




FIG. 2: A plot of the gluon distribution G(x) at Q 2 = 1.69 GeV 2 vs. Bjorken x. The solid curve is the published LO 
MSTW2008LO distribution and the (red) dots are the result of devolution from Q 2 = 5 GeV 2 to Q 2 = 1.69 GeV 2 . 

Appendix A: A Mathematica Laplace Inversion Algorithm 

Central to this numerical algorithm is the ability to write a closed form for the Pade approximant to 

1 



P(*) 



(Al) 



which was discussed in Eq. (33). 

Sidi [12., on p. 328, gives a closed form for the Pade approximant of iFi(l, f3, z). After some minor changes, we 
find that 



/'[ ^feM ,2JV- 1,2AM = 



(7) 



r(2N + j + /3 - 1) z 2Ar -J 



where 



(A2) 



r(/3+k) 



We now give the Mathematica algorithm that numerically implements Eq. (16) rapidly and accurately: 



NInverseLaplaceTransf ormBlockbeta[g_, s_, v_,betal_,twoN_, precision.] : =Module [ 
{Omega , Alpha , M , beta , p , den , r , num , hospital} , 

beta=Rationalize [betal , 0] ;prec=Max[precision,$MachinePrecision] ; 
M=2*Ceiling[twoN/2] ; Sp[pl_] :=Sum[(z~j)/Gamma[beta+j] , {j ,0,pl}] ; 

p=Expand[(Sum[(-l) "j Binomial [M, j] Gamma [j+M-l+beta] z~(M - j)*Sp[j-l], { j , 0, M}])]/ 



(A3) 
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Sum[(-l)~j Binomial [M,j] Gamma[j+M-l+beta] z"(M - j), {j,0,M}]; 
den=Denominator [p] ; r=Roots [den==0,z] ; Alpha=Table [r [ [i ,2] ] ,{i,l,M}] ; 
num=Numerator [p] ; hospital=(z~ (beta-1) num/D[den,z] ; 

Omega=SetPrecision [Table [hospital/. z->Alpha[ [i] ] ,{i,l,M,2}] ,prec+50] ; 
Alpha=SetPrecision [Table [Alpha [[i]] ,{i,l,M,2}] ,prec] ) ; 

SetPrecision[-(2/v) Sum [Re [Omega [ [i] ] g /.s -> Alpha[ [i] ] /v] , {i , 1 ,M/2}] ,prec] ] 



In the above algorithm, g = g(s), s = s, v — v, betal=/3, twoN=2iV in Eq. (16), and precision — desired precision 
of calculation. Typical values are twoN = 10 — 20 and precision ~ 30 — 100. The algorithm, which is quite fast, 
returns the numerical value of G(v). 

The algorithm first insures that M=twoN is an even number. It then constructs (in closed form) p, the Pade 
approximant of the function i-Fi(l, ft, z)/Y{ft), whose numerator is a polynomial in z of order twoN-1 and whose 
denominator is a polynomial in z of order twoN. It then finds r, the complex roots of the denominator, which are the 



a,, i.e., the poles of Eq. (16). Using L'Hospital's rule, it finds the residue Wj corresponding to the pole a*. At this 
point, all of the mathematics is symbolic. It next finds every other pair of (a,;, w,) to the desired numerical accuracy; 
they come consecutively, i.e., a% = a%, oj\ = u>2, = 014,, UJ3 = 0)4, etc. Finally, it takes the necessary sums, again 
to the desired numerical accuracy, but only over half of the interval i = 1,3,5,..., twoN, by taking only the real part 
and multiplying by 2. 

If g(s), the input to the algorithm, is an analytic relation and v and ft are pure numbers (from the point of view of 
Mathematica, 31/10 is a pure number, but 3.1 is not), then, for sufficiently high values of precision, you can achieve 
arbitrarily high accuracy. 

On the other hand, if g(s) is only known numerically, e.g., g(s) was obtained using numerical integration routines, 
the accuracy of inversion is limited by the need to only use relatively small values of 2N — in the neighborhood of 
2 — 8, limiting the overall accuracy to be in the neighborhood of 10~ 5 , which fortunately is ample for most numerical 
work. Typically, numerical integration routines are not accurate to better than ~ 10 -6 ; one can not use 10' a — which 
alternate in sign — that are larger than ~ 10 14 — 10 16 , which occur for relatively small values of 2N. Of course, this is 
not a limitation if g(s) is able to be expressed in closed form. 

We remind the reader that the algorithm can not be used for non-positive integral values of ft, because the exact 
inverse Laplace transforms are cither the Dirac S function or derivatives of it. 



Appendix B: Comparison with similar algorithms 

We had previously published two similar algorithms which we will designate as Algorithm I p] and Algorithm II 
[BJ. Algorithm I used the approximation 

1 /-Moo 2iV 



and thus is identical to the present work when ft — 1, as can be seen from Eq. ( 14 ). However, the algorithm for making 
the Pade approximant used in Ref. [I] is a factor of ~ 2.5 times slower than the algorithm that we give in Appendix 
|A") Although the majority of the computational time in both algorithms is spent in evaluating the real parts of the N 
complex g(cti/v), the new algorithm is slightly faster and is recommended. 
Algorithm II [BJ used the approximation 



G(v) « — / g(-) - 2 P(z 2 e z ,2N~l,2N), (B2) 

where P(z 2 e z , 2N — 1, 2N) is the Pade approximant of z 2 e z , whose numerator is a polynomial of order 2N — 1 and 
whose denominator is a polynomial of order 2N. It was specifically designed to make convergent original functions 
such as G(v) — 1/u 1- ^, ft < 0, and as such, required no knowledge of ft. However, as suggested earlier, an accurate 
numerical approximation to the value of ft is readily obtainable, with somewhat more effort. As an example, when 
using the original function G(v) = l/\/v of Ref. [BJ, the advantages of using our new algorithm is a speed factor 
of ~ 6, with a relative accuracy increase of O(10 -6 ), which is a significant gain; the disadvantage is the labor to 
determine ft. However, Algorithm II serves as a very useful numerical check that the g(s) indeed does asymptotically 
go as l/s? and that you are using an adequate numerical value for ft. 
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Appendix C: Hadamard Finite Part Integrals 



In this Section, we summarize some relevant mathematical details of (improper) Hadamard Finite Part ("parte 
finie") integrals, following the methods of Krommer and Ueberhuber [4], but modifying their notation somewhat after 
having adapted their work to our needs. We require these additional concepts when convolution integrals, such as 
Eq. ( |65[ ), involve kernels that are not regular functions, but rather are distributions. For the convenience of the reader, 
we again write the devolution relation of Eq. (651, i.e., 

G{v,Q 2 ) = K GG {w,T)G {v-w)dw+ K gf {w,t)F s0 (v - w) dw, (CI) 
Jo Jo 

recalling that negative r corresponds to devolution, i.e., evolving from Q\ to smaller Q 2 , since r was defined as 

Q 2 

r(Q 2 ,Q 2 ) = ~ I a s (Q ,2 )d\nQ' 2 , a s (Q 2 ) > 0. (C2) 

4n JQl 



The integral in Eq. (CI ) involving the kernel K G p{w, t) — > as r — > 0, so it is vanishingly small for all w and presents 
no problems for negative r. 

However, if r is negative and |r| <C 1, then the kernel K gg {w,t) must be a distribution in w which is almost a S 
function, implying that it can be written as 



K gg (w,t < 0) 



h(w) 



P<0, \0\ « 1, 



(C3) 



and h(w) is a smooth polynomial in w, with oc r. Clearly, whenever the overall exponent of w is greater than 1 in 
the denominator of Eq. (C3),when we insert it into the first integral on the r.h.s. of Eq. (CI), the integral diverges. It 
is this type of divergence problem, which occurs for all negative r in K gg (w,t), that we address in this Section. 



1. Giving meaning to the finite portion of a divergent definite integral 

To understand the concept of the Hadamard Finite Part integral, consider first the simple case of the divergent 
definite integral Jq(/9) for negative (3, defined as 

HP) = [^n dw > 

= lim / yr-— dw 



■^-}^ + wn- (C4) 



Equation (C4| shows that Jo can be broken up into two parts, the finite part f3~ 1 v 13 , which is called the Hadamard 
Finite Part integral, and a divergent part — lim,5^.o+ It is this finite part, with (i < 0, that is defined as 
the Hadamard Finite Part integral. We now introduce a new notation l/w~ l3+1 dw specifically for the Hadamard 
Finite Part integral, now writing 

W-£^^ = ^ /5<0- (C5) 
Next we generalize to a more complicated case, finding the Hadamard Finite Part integral of 

I f {fi) = f J^dw, /3<0, (C6) 
Jo w 

where f(w) is a general (Riemann integrable) function defined on [Q,v\. Let k = [— /3J, the floor of —fj, and define 
Tk(w), the Maclaurin polynomial expansion of degree k of f(w) as 

T k(w) = — w . (C7) 
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Since we will require that the Hadamard Finite Part integral defined in Eq. ( C6 ) inherit the Riemann integral properties 
of both linearity and additivity, we can rewrite Eq. ( C6 ) as 



Jo 



/H 



w 



dw 



f(w) - T k (w) 



1-/3+1 



1 



a J W-P+ 1 - 1 



dw. 



The first integral on the right-hand side of Eq. (C8|, 



I 



v f(w)-T k (w) 



dw, 



(C8) 



(C9) 



is readily seen to be an ordinary proper (or perhaps an improper) Riemann integral, so that we can now more simply 
write our definition of the finite part integral, Eq. (C8), as 



(0) = f 
Jo 



f_H 



dw = 



r f(w)-T k ( W ) a /«) r 1 , 

Jo v,-W + S^7o SP^*"' 



where the only Hadamard Finite Part integrals to be evaluated are of the form 

1 



£ 

■Jo 



-p+i-e 



dw, l = 0,l,...,k, fe=[-/3J, P < 0, 



(CIO) 



(Gil) 



which are readily evaluated using the results of Eq. (|C5j) . 

The short Mathematics program that follows, called intHFP[F,{w,0,v}] , evaluates the Hadamard Finite Part 
integral of Eq. (CIO I. For programming purposes, we have redefined the integrand F, using F = f(w)/w a+1 , whose 
variable w lies in [0, v], i.e., we have set a = —j3 > in the program. Note that one can alternatively use the form 
intHFP[F,{w,0,v},NIntegrate]to do numerical integration of the integral f Q (f(w) — T k (w))/w~^ +1 dw when v is 
input as a numerical quantity, if symbolic integration is neither possible nor desirable. 



Clear [intHFP] ; intHFP [F_, {x_, a_, b_} , int_: Integrate] := Module [{f , y, a] , sum , k , Tk} , 
{f ,a}=F/.f l_.x"al_.->{f 1,-al}; a = a-l; If [a<=0, Return [$Failed] ] ; 
k=If [IntegerQ [a] ,a-l,Floor[a]] ; 

If [a != 0, Return [Print ["Failed: Lower limit must be 0"]]]; 
sum=Total<aTable[((9 (Xii )f/(i!)/.x->0)*(-l/((a-l) b"(a-l),{i ) 0, k}] ; 
Tk=Total@Table[((<9 (;M jf/(i!)/.x->0) (x)~i, i, 0, k] ; 
int [(f-Tk)/x~(a+l) , {x, a,b}] +sum] /; 
(int===NIntegrate I I int==Integrate) ; Clear [x] 
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